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We present an efficient method for fast, complete, and accurate detection of unstable periodic orbits 
in chaotic systems. Our method consists of a new iterative scheme and an effective technique for 
selecting initial points. The iterative scheme is based on the semi-implicit Euler method, which has 
both fast and global convergence, and only a small number of initial points is sufficient to detect 
all unstable periodic orbits of a given period. The power of our method is illustrated by numerical 
examples of both two- and four-dimensional maps. 
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It has now been a widely accepted notion that unstable 
periodic orbits (UPOs) constitute the most fundamental 
building blocks of a chaotic system m . Theoretically, the 
infinite number of UPOs embedded in a chaotic invariant 
set provides a skeleton of the set, and many dynamical 
invariants of physical interest, such as the natural mea- 
sure, the spectra of Lyapunov exponents and fractal di- 
mensions, as well as other statistical averages of physical 
measurements, can be computed from the infinite set of 
UPOs in a fundamental way f^. In Hamiltonian systems, 
the quantum mechanical density of states in the semiclas- 
sical regime can be expressed explicitly in terms of UPOs 
of the corresponding classical dynamics ^J . The knowl- 
edge of UPOs is also of significant experimental interest 
because they provide a way to characterize and under- 
stand the chaotic dynamics of the underlying system [Q . 
All these call for efficient techniques for detecting UPOs 
in chaotic systems. 

Systematic detection of a complete set of UPOs of 
high periods embedded in a chaotic set even in situa- 
tions where the system's equations are known is, however, 
an extremely difficult problem. A fundamental reason is 
that the number of UPOs grows exponentially as the pe- 
riod increases at a rate given by the topological entropy 
of the chaotic set. The basic requirements for a good 
detection algorithm are, therefore, fast convergence and 
the ability to yield complete set of UPOs ||]. 

Recently, a general algorithm for detecting UPOs in 
chaotic systems was proposed by Schmclcher and Di- 
akonos (SD) Q who, for the first time, computed UPOs 
of high periods for systems such as the Ikeda-Hammel- 
Jones- Moloney map . The success of the SD method 
relies on a globally convergent iterative scheme: conver- 
gence to UPOs can be achieved, in principle, from any 
initial point. However, as we will discuss shortly, this 
method is not very efficient from the standpoint of con- 
vergence, neither does it provide a satisfactory test for 
the completeness of the detected UPOs. As a matter 
of fact, for the Ikeda-Hammel-Jones-Moloney map, only 
UPOs of periods up to 13 were reported in Ref. 1^, and 
one of the UPOs of period 10 was not detected. 

The aim of this Letter is to present an efficient method 
for detecting UPOs in general chaotic systems. Our 
new iterative scheme is based on the semi-implicit Eu- 



ler method S and has the following favorable proper- 
ties: near an orbit point it exhibits a fast convergence 
similar to that of the traditional Newton-Raphson (NR) 
method, while away from the orbit points it is similar 
to the SD method and, therefore, is globally convergent. 
Another key ingredient of our method is that we select 
initial points based on the observation that using orbit 
points of UPOs of other periods to initialize the search for 
UPOs of a given period is much more effective than using 
randomly selected points in the phase space or in the at- 
tractor. We find, in most cases, it is sufficient to use only 
orbit points of period p — 1 in order to detect all UPOs of 
period p. With such a strategy, we are able to compute 
UPOs for, say, the Ikeda-Hammel-Jones-Moloney map, 
of periods up to 22 for a total of over 10® orbit points 
using roughly the same amount of computation required 
by the SD method to compute all UPOs of periods up 
to 13 that have less than 6000 orbit points |^. Due to 
its efficiency, our method allows us to compute UPOs in 
higher-dimensional systems, which we illustrate using a 
four-dimensional chaotic map. 

We begin by describing the NR and the SD methods. 
Consider an iV-dimensional chaotic map: x„+i = f(x„). 
The orbit points of period p can be detected as zeros of 
the following function: 



g(x) = f(P)(x)-x, 



(1) 



where f*^P)(x) is the p-times iterated map of f(x). The 
process of finding zeros of g(x) usually begins with the 
choice of initial point xq followed by the computation of 
successive corrections: Xnow = Xom + 6x, which converge 
to the desired solution. In the NR method, the correc- 
tions are calculated from a set of N linear equations: 



J(x)5x = g(x) 



(2) 



where J(x) — dg/dx is the Jacobian matrix. The NR 
method has excellent convergence properties, approxi- 
mately doubling the number of significant digits upon 
every iteration, provided that the initial point is within 
the linear neighborhood of the solution. While it is rel- 
atively easy to find suitable initial points for very small 
periods (using, for example, uniform grid, iterations of 
the map, or random number generator), the method be- 
comes impractical for UPOs of high periods because the 



volume of the basin from which Xq can be chosen de- 
creases exponentiahy as the period increases. On the 
other hand, in the SD method, the corrections are deter- 
mined as foUows: 



<5x = ACg(x) , 



(3) 



where A is a small positive number and C is an iV x iV 
matrix with elements Cij G {0, ±1} such that each row or 
column contains only one element that is different from 
zero. With an appropriate choice of C and a sufficiently 
small value of A the above procedure can find any pe- 
riodic point of a chaotic system. The main advantage 
of the SD method is that the basin of attraction of each 
UPO extends far beyond its linear neighborhood, so most 
initial points converge to a UPO. In fact, the basins of at- 
traction of individual orbit points completely fill a region 
in the phase space, and any initial point in this region 
converges to an orbit point. 

Schmelcher and Diakonos tested their method by com- 
puting the UPOs for the Henon map and other simple 
maps, for which the UPOs are known from methods spe- 
cific to these maps [||. They also applied the method 
to the Ikeda-Hammel-Jones-Moloney map, for which no 
special technique for computing UPOs was previously 
available. The method appears to be particularly useful 
when detecting for each period the least unstable peri- 
odic orbits iQ. However, if the goal is to determine com- 
plete sets of UPOs of increasingly higher periods, the SD 
method becomes inefficient due to the following two rea- 
sons: (i) the convergence rate of Eq. (ph is much slower 
than that of the NR method, so it takes significantly 
more steps to reach the desired accuracy [|ll[; and (ii) 
even though the SD scheme is globally convergent, the 
basins of attraction of individual UPOs are interwoven 
in a complicated manner, so it is extremely difficult to 
determine which initial point converges to a particular 
UPO. Because of this difficulty, the SD method cannot 
guarantee the detection of all UPOs of a given period. 

To overcome the problem of slow convergence, while 
retaining the global convergence property, we propose the 
following iteration scheme: 



[l/3<7(x) - CJ(x)]<5x = Cg(x) 



(4) 



where ^(x) = ||g(x)|| > is the length of the vector, 
/3 > is an adjustable parameter, and C is the same 
matrix as in Eq (Q). In the vicinity of an UPO, the func- 
tion (7(x) tends to zero and the NR method is restored. In 
fact, it is straightforward to verify that the above scheme 
retains the quadratic convergence. On the other hand, 
away from the solution and for sufficiently large values of 
/3, our scheme is similar to Eq. (||) and thus almost com- 
pletely preserves the global convergence property of the 
SD method. This similarity is easily understood, since 
Eq. (13) is the Euler method with step size A for solving 
the following system of ODEs: 



s^-^«- 



(5) 




On the other hand, Eq. (0) is the semi-implicit Euler 
method |^ with step size h = 1/ j3g{'K) for solving the 



FIG. 1. Shown with thick lines are the basins of con- 
vergence of (a) the Newton-Raphson (NR) method, (b) the 
Schmelcher-Diakonos (SD) method with < A < 0.3568 and 
C = 1, and (c) our method with /3 = 4.0 and C = 1 to the 
zeros of a function cos(x^) in the interval (—3,3). Arrows 
indicate the direction of convergence and large dots are the 
zeros to which the methods converge. 



same system of ODEs. Consequently, with sufficiently 
small step size, both methods closely follow the solutions 
to Eq. (P) and thus share the global convergence property. 

To illustrate and contrast the convergence properties 
of the NR, the SD, and our methods, we consider the 
following simple example: finding zeros of the function 
g{x) = cos(a;^) in the interval (—3, 3). The basins of con- 
vergence for each method are shown in Fig. n^ with the 
thick arrows. The NR method converges to all six zeros 
and the basins are essentially within the linear neigh- 
borhood of each point [|2|. The SD method converges 
to the solution g{xo) = if < A < 2/5' (a;o) and 
C = —sig\Y{g' {xq)) . Diagram (b) in Fig. [^ shows the 
basin of convergence for < A < 0.3568 and C = 1. 
Obvious is the global character of convergence to zeros 
with negative function derivatives, while zeros with posi- 
tive derivatives serve as basin boundaries. With C = — 1 
the convergence directions are reversed. The result of 
applying our iteration scheme with (3 — 4.0 and C = 1 
to the same function is shown in the diagram (c). We 
see that, as in the NR method, all zeros have basins of 
convergence. However, the basins of zeros with negative 
function derivatives cover most of the interval, while the 
basins of other zeros, as well as the intervals between 
basins, are reduced and become smaller with increasing 
value of (3. Therefore, our scheme combines the efficiency 
of the NR method with the global character of the SD 
algorithm. 

Another important ingredient of our method lies in the 
selection of initial points: we find that the most efficient 
strategy for detecting UPOs of period p is to use UPOs 
of other periods as initial points. This is understandable, 
since orbit points cover the attractor in a systematic man- 
ner, which reflects the foliation of the function f ^^^ (x) and 
its iterates. In cases of the Henon and the Ikeda-Hammel- 
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FIG. 2. Number of detected orbits for different periods in 
the Ikeda-Hammel-Jones-Moloney attractor given in Eq. (pi). 
Solid dots indicate the values of (3 above which we always 
detect a maximum number of UPOs for each period. 

Jones-Moloney maps, we are able to detect all UPOs of 
period p using only orbit points of period p — 1, provided 
that period p — 1 orbits exist. In more complicated cases 
of higher-dimensional maps, this simple strategy leaves a 
small fraction of UPOs undetected [Ol. However, in all 
cases, we are able to find these UPOs using period p+1 
points (first we use incomplete set of period p orbits to 
find period p+1 points and then use them to complete 
the detection of period p orbits) . The main advantage of 
using orbit points of neighboring periods as initial points 
is that once we establish the strategy for smaller periods, 
it works in a similar manner for the detection of UPOs 
of large period. This allows us to claim with confidence 
that we detect all UPOs of increasingly longer periods 
for general multi-dimensional chaotic maps. 

We now apply our method to detecting UPOs for the 
following Ikeda-Hammel-Jones-Moloney map |0] : 



x' = a + b{x cos <j) — y sin ( 
y' — b{x sin (f> + y cos (f)) , 



(6) 



where </) = k—ij/^l+x'^+y'^), and the parameters are cho- 
sen such that the map has a chaotic attractor: a — 1.0, 
h = 0.9, k = 0.4, and -q = 6.0. Detection of UPOs pro- 
ceeds as follows: UPOs of period 1 and 2 are quickly 
found using several initial points on the attractor. Start- 
ing from p = 3 we use only orbit points of period p—1 as 
initial points. We choose C from the set of five matrices 
{Cfcjfc = 1,...,5} provided in Ref. [|j, where Ci = 1 
is the identity matrix. The iteration sequence computed 
from Eq. (0) is terminated when it either converges to an 
orbit point or leaves the chaotic attractor. The average 
number of iterations increases linearly with /3, which is 
understandable since ||5x|| Ki 1//3 for large (i and away 
from an orbit point. However, a small fraction of initial 
points produces very long sequences which neither con- 
verge to an UPO nor leave the attractor. In order to 
limit the amount of unproductive computation, we set 
the maximum number of iterations to 4-6 times /3, which 



TABLE I. Number of distinct UPOs, n(p), and the total 
number of orbit points or period p, N(p), for the Ikeda at- 
tractor given by Eq. (pi). Note that N{p) also includes orbit 
points whose periods are factors of p. 



P 



<P) 



N(p) 



14 
15 
16 
17 
18 
19 
20 
21 
22 



317 

566 

950 

1646 

2 799 

4884 

8404 

14 700 

25 550 



4511 

8517 

15 327 

27 983 

50 667 

92 797 

168 575 

308 777 

562 939 



is sufficient for the majority of iterates to be terminated 
properly. The quadratic convergence of our scheme al- 
lows us to achieve, without much computational effort, 
accuracy limited only by the computer round-off error. 
Once the sequence converges to an orbit point, we check 
whether it belongs to a yet undetected UPO, and if so, 
we compute the rest of the orbit points by iterating the 
map and refining the solutions with a couple of NR steps 
[we simply set (3 — in Eq. (^) ] . 

Figure shows the number of detected UPOs of peri- 
ods 10 through 18 using different values of /? in the range 
from 10 to 3000. Note that for every period there exists 
a value /3 — Pmin{p) above which we are guaranteed to 
find a maximum number of UPOs. This feature of our 
scheme strongly suggests that the detected orbits con- 
stitute a complete set of UPOs for each period. Since 
/3min(p) is approximately proportional to e"P, where a is 
a positive constant, we can estimate the value of /? neces- 
sary to find all UPOs of increasingly longer periods. The 
numbers of the UPOs for periods up to 13 agree with 
those of Schmelcher and Diakonos g except for period 
10, where we have detected one additional orbit. The 
number of orbits of periods 14 through 22, which were 
not reported previously, are given in Table B. 

If we monitor the number of orbits detected with differ- 
ent matrices C, we note that, for a wide range of values 
of /3, after we use identity matrix Ci, only a few UPOs 
remain undetected. For example, with /? — 5000 and 
C = Ci in Eq. (Q), our method detects 14 699 orbits 
of period 21, and only one new orbit is detected with 
C = C2. To understand this feature of our method, 
which is common to all the maps tested, we show in 
Fig. y, for the chaotic attractor in Eq. (H), the num- 
ber of period 13 orbits detected with Ci (solid dots) 
and the number of additional orbits detected with Cfe, 
k = 2,..., 5, (triangles). For 100 < /3 < 1000, almost 
all UPOs are detected with C being the identity matrix. 
At larger values of (3 the number of thus detected orbits 
decreases, but the remaining orbits are always detected 
with other matrices. For /3 > 10^ the numbers converge 
to those of the SD iteration scheme, where about half of 
the orbits are detected with Ci and the other half with 
C2 and C3. This behavior of our scheme follows directly 
from the convergence considerations of Fig. |l| and results 
in a greatly improved efficiency compared to either the 




FIG. 3. Detection of UPOs of period 13 in the 
Ikeda-Hammel-Jones-Moloney attractor. The number of or- 
bits detected with Ci is shown with sohd dots, while triangles 
represent the number of additional orbits detected with C^, 
k = 2, . . . , 5. The total number of detected orbits is shown 
with open circles. 



NR or the SD methods. 

Finally, we briefly describe the performance of our 
method for other maps. In case of the Henon map 
our algorithm works extremely well and, for the stan- 
dard parameter values of (a, 6) = (1.4,0.3), detects all 
UPOs up to period 29 with /3 < 500, C = Ci and 
C2, and using for initialization only orbit points of pe- 
riod p — 1. We have also applied our algorithm to de- 
tecting UPOs in the following four-dimensional map: 
Two coupled Ikeda maps with coupling in the form: 
0(1,2) = ^-'7/(1+2^^1,2) +yfi,2)) + 27re(a;(2,i)-a;(i,2)), and 
the parameters are chosen such that the system has two 
positive Lyapunov exponents. We estimate the topolog- 
ical entropy in this system to be /it ~ 1.6, and thus the 
number of orbits grows extremely fast with increasing or- 
bit length. We have detected complete sets of UPOs up 
to period 7 with /3 < 1000. We have found that the relia- 
bility of the algorithm was not affected by the increased 
dimensionality of the system. Even though the number 
of possible matrices C in four dimensions is 384, only a 
dozen of them are needed to detect all UPOs. The neces- 
sary set of matrices C can be selected empirically when 
detecting short UPOs and then used in the detection of 
longer orbits. 

In conclusion, we have proposed an efhcient algorithm 
for the detection of UPOs in chaotic systems and have 
successfully detected large number of UPOs in several 
two- and higher-dimensional maps. Our method allows 
for a verification of the completeness of the detected or- 
bits and high accuracy limited only by the round-off er- 
ror. 
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